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We develop a semiclassical density functional theory in the context of quantum dots. 
Coulomb blockade conductance oscillations have been measured in several experiments using 
nanostructured quantum dots. The statistical properties of these oscillations remain puzzling, 
however, particularly the statistics of spacings between conductance peaks. To explore the role 
that residual interactions may play in the spacing statistics, we consider many-body systems 
which include electron-electron interactions through an explicit density functional. First, we 
develop an approximate series expansion for obtaining the ground state using the idea of the 
Strutinsky shell correction method. Next, we relate the second-order semiclassical corrections 
to the screened Coulomb potential. Finally, we investigate the validity of the approximation 
method by numerical calculation of a one-dimensional model system, and show the relative 
magnitudes of the successive terms as a function of particle number. 

PACS numbers: 73.23.Hk, 05.45.Mt, 71.10.Ca, 71.15.Mb 

I. INTRODUCTION 

A recurring problem in modern physics is how to add quantization effects to a basically successful macro- 
scopic theory. This question arises particularly in the semiclassical regime — large quantum number — where 
the quantum effects are often corrections to the essentially classical macroscopic physics. Perhaps the best 
known example starts with the Thomas-Fermi theory of the atom |^], which is macroscopic in essence, and 
then evaluates the contribution of electronic shell structure to the ground state energy . The very natural 
result that the shell contribution is given by the quantized levels of the self-consistent Thomas-Fermi potential 
has been used extensively [^|-^. However, it has only been in recent decades, starting with the work of V. M. 
Strutinsky [^-^ , that a systematic way of answering the recurring general problem has been developed. 

Our own immediate interest is in quantum dots — small electrically conducting regions in which the quantum 
properties of the confined electrons are important pO[ | — and our aim here is to treat quantum corrections to the 
ground state energy of these dots by further developing the Strutinsky method. 

Quantum dots can be formed, for instance, by gate depletion of a two-dimensional electron gas (2DEG) in a 
GaAs-AlGaAs heterostructure. Because of the high quality of this material and interface, the mean free path 
of the electrons far exceeds the size of the quantum dot. One can view an electron as propagating ballistically 
within a confining potential created by electrostatic gates patterned on the surface of the heterostructure. For 
transport measurements, dots can be coupled weakly to leads; when the conductance of each lead falls below 
2e^//i, electron transport through the dot occurs only by tunneling, and the number of particles within the dot 
becomes quantized. In this regime, the conductance is suppressed due to the electrostatic energy associated with 
a localized charge an effect known as the Coulomb blockade (first reported in 1951 by Gorter |l^). When 
the dot potential is tuned by a gate voltage so that adding one electron costs no energy, a large conductance 
peak appears [|l^ — though the electrons must still tunnel from the leads, there is no additional electrostatic 
barrier to conduction. Sweeping such a gate voltage produces periodic Coulomb blockade oscillations of the 
conductance through dots |p^ , |l2| , |l3| . 

The Coulomb blockade is a classical effect observable in a broad temperature range, ksT < e^/C, where C is 
the total dot capacitance. Over most of this range, both the spacing and height of the peaks is constant — the 
spacing is e^/C and the height is given by the resistance of the two tunneling barriers acting in series. However, 
there is a low temperature regime below a few hundred millikelvin for which ksT < A, where A is the mean 
single-particle level spacing of the isolated dot. There, quantum interference and coherence become important 
|fiO|. The Coulomb blockade peaks grow as the temperature decreases, and novel fiuctuation properties emerge 
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involving both the peak heights ]T^~|T6t and spacings [^7|-|T] . The spacings give information about the ground 
state energies while the heights involve the magnitude of the wavefunction near the levels. 

With some success, random-matrix-theory (RMT) based approaches have been used in order first to 
predict psf and then to explain the statistical properties found. In the simplest approximation [T^ ], known as 
the constant interaction model, the ground state energy of the dot is expressed as 

i=l 

where N is the number of electrons in the dot and are single-particle energies. The first term is the classical 
charging energy; the second is the total energy of a system of non-interacting quasi-particles. Supposing 
that the single-particle classical dynamics within the dot is chaotic, the Bohigas-Giannoni-Schmit conjecture 
applies ||2^ and implies that the single-particle quantum properties follow RMT. It is well-established that 
RMT predicts repulsion amongst the Ci and Gaussian random behavior in the eigenfunctions psl . Within this 
model, the conductance peak height statistics are in good agreement with experimental results [|l4| , p^ after also 
incorporating non-zero temperature effects and interference modulations due to periodic paths coupled to the 
leads p3| , p6|j27| . On the other hand, it is found that the fluctuations in the peak spacings are considerably 

larger than the predictions ||l^,^,0 , and there is no evidence for the level repulsion or electron spin degeneracy 
expected from a single-particle-like approach | |30|] . These discrepancies between the predictions of the constant 
interaction model and the observations point to the need for a quantum treatment of the electron-electron 
interactions, and, in particular, have triggered a number of studies based on Hartree-Fock calculations 31-33|, 



or density functional theory in the local density approximation (LDA) [p4 35 1. 

In nuclear physics, it has been recognized for some time that the dependence of nuclear ground state prop- 
erties on the particle number can be viewed as the sum of two contributions, a "smooth part" which varies 
smoothly with the particle number and an "oscillatory behavior" producing the shell structure referred to as 
shell corrections. A similar decomposition is possible for any finite-size interacting fermion system. The smooth 
part comes basically from the bulk energy per unit volume integrated over the finite-size system, and the oscil- 
lating contributions come from quantum interference effects explicitly caused by the confinement. By supposing 
that the smooth part is known while the unknown oscillatory contribution is a correction, Strutinsky proposed 
in the 60's a physically motivated approach to, and an efficient way of calculating, the shell corrections 

Strutinsky's shell correction method is essentially a semiclassical approximation. It rests on the fact that the 
number of particles in the system considered is large, rather than on the interaction between the particles being 
weak. (One must, of course, work in a regime where the smooth starting point is basically valid.) Since the 
quantum dots in which we are interested contain on the order of a hundred electrons, they are a perfect place to 
apply the Strutinsky method. However, before doing so for a particular, realistic, two-dimensional geometry, we 
shall in this paper limit ourselves to a formal discussion of this method in conjunction with a one- dimensional 
illustrative example. In spite of the literature existing on this subject we find it useful for two main purposes. 
First, the discussion and resulting expressions are noticeably simpler for quantum dots than for nuclei. This 
occurs because the existence of a smooth confining potential in the dot means that gradient corrections to the 
smooth density are not needed to confine the system at the zeroth-order (classical-like) approximation. The 
effect of these gradient terms can therefore be included in the first-order "shell" corrections, simplifying both 
the zeroth-order calculations (no gradient terms) and the first-order terms (no corrections to the Weyl part). 

Our main purpose, however, is to take advantage of the fact that we use density functional theory rather 
than Hartree-Fock as a starting point, the former being presumably better suited to deal with the long-range 
Coulomb interaction present in quantum dots than the latter. This allows us to discuss in detail the second-order 
"residual interaction" terms of the Strutinsky method. By residual interaction we mean the weak interaction 
between Landau quasi-particles that comes from dressing the bare electron added to the quantum dot. In 
particular we will show, and illustrate, how these terms are related to the screened Coulomb interaction. 

The remainder of the paper is organized as follows. The Thomas- Fermi and density functional theories are 



summarized in the next section, establishing our notation. Section III contains the Strutinsky method applied 
to density functional theory. This is the core of the paper; in particular, the relation of the second order terms 
to the screened Coulomb potential is derived. Section ^ recalls how the residual interaction terms contribute 
to conductance peak spacing distribution. Section ^ compares the whole approximation scheme to numerical 
calculations of a simplified model: interacting electrons in a one-dimensional quartic oscillator. Finally, we 
comment on the relationship between the Strutinsky development and the constant interaction model, and 
possible applications of the method. 
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II. DENSITY FUNCTIONAL THEORY 



The Hohenberg-Kohn theorem states that for a system of interacting electrons in an external potential, 
V;,xt(r), there exists a functional, .7-"hks["]j of the density of electrons, n{r), such that: i) the density, ng{r), 
corresponding to the ground state of N particles is an extremum of ^hks under the constraint that the total 
number of particles, 

N[ng] = / drngir) , (1) 



is fixed, and ii) J-HKs["-g] is the total energy of the system. The explicit form of the Hohenberg-Kohn-Sham 
functional is not known p6| , |37[ |. In practice, one must be satisfied with approximations. We describe here first 
a generalized Thomas-Fermi approach and, second, the case when an explicit form of the density functional is 
assumed. 

A. Generalized Thomas- Fermi Approximation ([Semi] Classical Level) 

It is convenient to view the density functional as the sum of three parts: a classical charge contribution, the 
kinetic energy, and the unknown exchange-correlation functional which accounts for the balance js^] . The first 
part is simple: the energy of a system of classical charges confined by an external potential, Vcxt, is 

£[n] = £cxt[n] + fcouib^] (2) 

where 



£cxt[n] = J n(r)Kxt(r)dr 

^cou.[n] = ^ / / ^P^drdr' . (3) 



For the kinetic energy, in the Thomas- Fermi approach the Pauli exclusion principle is introduced semiclassi- 
cally by employing the idea that one quantum state occupies a volume (27r?i)'' in phase space. This implies that 
if many electrons want to be at the same place, they can do so only by increasing their kinetic energy. This 
gives the Thomas- Fermi approximation to the kinetic energy part of the density functional, 7tf["-]j expressed 
as 

1 f f p2 



^ ' (2TrhY J \ 2m/ ^ 



tTp{n) ~ / e{i')di' 
Jo 

TtfM = J tTF{n{r))dr (4) 

where d is the dimensionality of the system, Q is the Heaviside step function, ixF is the kinetic energy density, 
and z^(e) is the number of states per unit volume with energy less than e. A factor 2 in i>{e) is required if the 
electron spin degeneracy is taken into account. 

Finally, the effect of exchange and correlation is included through a term £xc['t-]- In practice, an explicit form 
for this functional must be taken. For example, if the electron density is a sufficiently slowly varying function 
of position, one can approximate fxcl?^] by taking the exact results for the uniform electron gas at the local 
density, the well-known local density approximation (LDA). 

Within this approximation, then, the density functional is 

^GTF [n] = Ttf [n] + Stot [n] (5) 

where 

£tot[n] = fcxtM + fcouiN + £Kc[n\ ■ (6) 
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The ground state energy and its electron distribution are obtained by minimizing JFqtf under the constraint 
(Qj), yielding the self-consistency equation 



Sn 

with the effective potential 



[noTF] (r) + Vcs [noTp] (r) = Mgtf (7) 



KffM(r)^^N(r). (8) 



Notice that to make use of Eq. (|^), one must have an explicit form for £xc in order to take the functional 
derivative in Eq. (||). 

We call this approach "generalized-Thomas-Fermi" (GTF) because it uses the Thomas-Fermi approximation 
for the kinetic energy but retains macroscopic aspects of exchange and correlation. In particular, the short-range 
effects of exchange-correlation in a uniform system can be included. 



B. Kohn-Sham Equations (Quantum Mechanical Level) 

In standard implementations of density functional theory (DFT), the kinetic energy is treated quantum 
mechanically rather than (semi)classically as in GTF. To accomplish this, one introduces N orthogonal functions 
{01 (r), (^jv(r)} such that the electron density is defined as 

N 

n(r)^^|0,(r)p (9) 



and the kinetic energy is 



N 



N 

. y ^|V0.(r)pdr. (10) 
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Thus, the density functional becomes 

^DFT W = TdFT W + £tot[n] ■ (11) 

To find the ground state energy, one minimizes J-uft with respect to the functions (f>i (r) under the constraints 



yi0,(r)pdr = l, z = l,...,iV. (12) 



The result is a Schrodinger equation 

— V2 + y,ff[n](r)j 0,(r) (r) , i = l,...,N (13) 

where ei, ■■■,eN are Lagrange multipliers and the effective potential is again defined by Eq. (|^); these are the 
Kohn-Sham equations |^^. The Eqs. (^), (||), and (|l3|) are the set of self-consistent equations for finding the 
electron density, tt-dft (r) , and then the ground state energy E'dft = -^I^ft ['^dft] • 

As in the discussion of the GTF above, in order to actually solve the Kohn-Sham equations, an explicit form 
for the exchange-correlation functional is required. The simplest case is when 5xc is an integral over space of 
a function (not functional) of the local density; this is the well-known local-density approximation (LDA). But 
other more complicated explicit forms are possible, for example the generalized gradient corrections to LDA p8|. 
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III. APPROXIMATE GROUND STATE ENERGY: 
THE STRUTINSKY ENERGY-CORRECTION METHOD 



A. Expansion of the DFT Ground State Energy 



In this section we develop an approximation to i?DFT [jt^dft] starting from the solution of the generalized- 
Thomas- Fermi equation, ncTF- The main motivation is to develop a physical interpretation of the difference 
between these two approaches in finding the ground state energy. In addition, the approximation is of interest 
numerically for large problems since it involves a self-consistent solution of only the GTF equation rather than 
the more involved Kohn-Sham equations. We use the method introduced by V. M. Strutinsky originally 
in the context of a Hartree-Fock rather than density-functional approach. His method describes the interacting 
system self consistently, first with the quantum interference effects turned off, and then by introducing them 
perturbatively. As discussed in the introduction, the idea is to add the "oscillatory" effects caused by interference 
in the confined system to a "smooth" essentially macroscopic description — these effects are essentially the Friedel 
oscillations [ ^9| familiar in the context of impurities or surfaces. 

To study the role of quantum interference effects in the DFT ground state energy, we will first show that the 
generalized- Thomas-Fermi result is a [semi] classical approximation to the DFT energy. The GTF approximation 
does, of course, contain some quantum mechanics — notably the Pauli exclusion principle which gives rise to 
the Fermi surface — and so is not truly classical. But only the simplest local quantum effects are present in 
GTF rather than the effects of interfering paths that one expects in a true semiclassical theory, hence our 
characterization of GTF as " [semi] classical" . 

To see this clearly, we introduce a convenient notation adapted from the semiclassical treatment of single- 
particle problems: it is customary there to express the density of states as a sum of a smooth term slowly 
varying in energy, called the Weyl part, and a term which varies rapidly in energy (on the scale of the mean level 
separation), called the oscillatory part [^]. For a system governed by the Hamiltonian H[V] = Tp^ /2m + V{y), 
where the potential is as yet unspecified, one can define the probability density of N independent particles 



N 

n[V^](r)=5]|^,(r)l2 (14) 

4=1 

in terms of the eigenstates {(t)i] of H . We also define the Weyl part of n\V] by 

n'^mir) EE 1 e [m^ - pV2m - V{r)] dp (15) 

where fj}^ must be chosen so that N = J (r)dr. Note that n'^[y](r) is smooth in that it neglects quantum 
fluctuations in much the same way that the GTF approximation does. With this notation, one can derive the 
useful relation 



6n 

Indeed, using 



[n^[V]] (r) +Vir) = fi^ . (16) 



(57tf 



n](r)-e(n(r)) , (17) 



Eq. (|T^) reads e(n'^[y](r)) = fi^ — V{r). Applying the function ly introduced in Eq. (||) to both sides of the 
equality gives the definition of [V] (r) , Eq. ( 1^) . 

Recalling that T4ff[n] is defined as the variational derivative of £tot (to be completely clear, it is not the inverse 
of n[Kff]), we see that the self-consistency equation (^ which defines uuft is 

'^DFT(r) = '^[Vcff[nDFT]](r) • (18) 

Similarly, Eq. (|^) which defines ncTF can, in applying the above prescription, be put in the form 

nGTF(r) = [Ks [t^ctf]] (r) . (19) 
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These equations do not signify that ncTF is the Weyl part of Udft; however, they do indicate that if one 
neglects the quantum interference terms (i.e. the difference between the exact particle density and its Weyl 
part), then the definitions of uqtf and tt-dft become equivalent. It is in this sense that tigtf is the [semi] classical 
approximation of noFT- 

Supposing nGTF(r) and -Eqtf = •^GTFi'^GTF] known, we now seek to evaluate the corrections to the Thomas- 
Fermi energy, 

AE = EuFT — Eqtf , (20) 

up to second order in 

Sn = riDFT - ncTF ■ (21) 
For this purpose, we first introduce the quantities 

n(r) = n[l/cff[nGTF]] (r) 
n^(r) = [V;ff[nGTF]] (r) = '^GTF(r) 

n°"'{r)^fi{r)~fi^{r) . (22) 

Note that once 14ff [tigtf] is known, all of these can be computed through the diagonalization of the known 
single particle GTF Hamiltonian. As is well-known, the sum of the eigenvalues of the Kohn-Sham equations, 
£ip[]/] = J2i does not give the total energy of the N particles because of double counting of the interaction 
energy, but rather 

^^DFT = flp[K!ff["-DFT]] - y"drVoff [nDFT](l")"-DFT(l") + '^^tot ["-DFt] ■ (23) 

To proceed further, we use the relation proved in the appendix 

£ip[V + SV] - £ip[V] ~ i y" SV{r){n{r) + n'{r)) dr (24) 

where n(r) = n[y](r) and n'(r) = n[V + SV](r) and which is correct through second order in the changes. 
Upon inserting V = Fcff [hgtf] and V + SV = T4ff[nDFT], and thus n = h and n' = tidft, the first term on the 
right-hand-side of Eq. ([2^) becomes 

fip [Kff[nDFT]] = Sip [Kfr[^^GTF]] + dr SV,s{r){nj,FT{r) + n(r)) (25) 

where SVcs = Vcff[nDFT] — K!ff["-GTF]- Similarly, the second term in £'dft['^dft] is 

J VcS [udft] (r)'^DFT {r)dr = J T4ff [f^GTF] (r) "^gtf (r) dr 

+ J 14ff[nGTF](r)(5n(r)rfr (26) 

+ j (^Kff (r) ?^DFT(r) dr . 

Finally, the third term is 

Aot[nDFT] = ftot[nGTF] + J (Kff [^gtf] (r) + SV,s{r)/2)Sn{r)dr (27) 

with corrections which are third order in Sn. 
Combining all the terms together, we obtain 

EuFT — flp[Voft[nGTF]] - y nGTF(r) Voff [nGTF](r) dr + £tot[nGTF] 

+ ^ /^Feft(r)n°-(r)dr . (28) 
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In order to express directly the difference between the DFT and GTF ground states, it is convenient to use 

Egtf [ncTp] = Ttf [^ctf] + ^ tot [«gtf] (29) 

for ftot- In order to simphfy the last term in Eq. (p8|), note that n"^"^ is of order Sn, and that therefore one only 
needs the first order variation of the effective potential, SVcs{r) — J {SVes/6n)[nGTF]{Tc,r')Sn{r')dr' , to obtain 
-E'dft correct through second order. Thus the final expression for the Strutinsky energy correction is 

£;dft ^ Egtf + AE^^^ + AE^^^ (30) 

where the first and second order correction terms are 

AE^'^^ = fip[Foff[nGTF]] - J nGTF{r)Vcs[nGTF]{r)dr ~TTF[nGTF] (31) 

^^^'^ ^IJ I ^"'^'■^ ^t''^^'']^'"' ""'^ '^''^'■'^ ■ ^^^^ 

In this approach, the DFT ground state energy is, then, the sum of a classical contribution — the generalized- 
Thomas- Fermi result -Egtf — and two quantum contributions — AE^^^ and Ai?*^^^. We now discuss and interpret 
these two correction terms. 

B. Interpretation of the First-Order Corrections 

The first-order correction is simply the oscillatory part of the single particle energy for a system of N electrons 
evolving in the potential Vi3ff[nGTF]- Indeed, the Weyl part of fip[V^] is 

^ZiV] - tAu. [(^ + V{r)) eU""' - 1^ - V{r))dpdr (33) 



^P^^~ i2TTh)d J \2m ^ 'J V 2m 

where /z^ is fixed by iV = J n}^ [V]{r)dr. Separately integrating the kinetic and potential energy terms for 
V = Kjff[nGTF], one obtains 

r /•m"' -Kff[nGTF](r) J p 

£^ [Feff [ncTF]] = J drj e—de + jn"^ [Kff Ktf]] (r) T/eff [hctf] (r) dv . (34) 



In the first term one recognizes the Thomas- Fermi kinetic energy, '7tf['^gtf], while in the second term 
[Vi3ff[nGTF]] = nGTF- Thus the first-order Strutinsky correction is 



w 



AS(i) = £ip [yeff [ncTF]] - £Z [Kff[nGTF]] = f [Ks^tf]] ■ (35) 

The leading quantum corrections to GTF are found, then, by quantizing the single-particle levels in the GTF 
self-consistent potential: this is a very natural result which, in fact, was used extensively in atomic and nuclear 
physics ||-^ before it was first justified by Strutinsky 0,^. 

C. Interpretation of the Second-Order Corrections 



The second-order correction, Eq. (|3^), requires further work: this form is not useful because it expresses 
Ai^^^' as a function of the unknown 5n. A second equation is necessary for us to determine 5n. Note that this 
is not the case for AE^^^ which is written completely in terms of tigtf- 

The required second equation is obtained by relating 5n to the oscillatory part of n = n[yoff [tjctf]] which, 
of course, is known since it depends only on ugif- We start with the two equations 

(5'7xF 

-^[^GTf] + ^^({[nGTF] = MGTF (36) 

-T^["dft] + Kff[^^DFT] = Mdft • (37) 
on 
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The first equation here is the definition of ncTF, and the second one follows directly from the general relation 
(p^. Now expand tt-dft about ncTF in the second equation and subtract the first one from it. In the term 
involving Kff, 5n appears. However, in the kinetic energy term, the density difference is n^^rp — tigtf = 
("-DFT — ^^ctf) — (JT-DFT — 't-dft) = ^''^ " ^^DFT- ^0 close the equation we must relate rippx to ncTF- This 
is possible because in an equation for dn, which is by definition first order in corrections, only the first order 
part of the other quantities need be kept. Thus, we can approximate n^prp by similarly expanding tidft about 
rtGTF, yielding 

"DFT = ("DFT - ngpT) - " ^gtf) = n"'" ■ (38) 



The combination of these results gives the closure equation 



J dr^(r,r')<5n(r') + J dr^^(r,r') {Snir') - h°'%r')) = 



(39) 



where A/^ = /if^^p — Mgtf is fixed by the condition J 6n{r)dr = 0. This is an integral equation for Sn in terms 
of GTF quantities. If a numerical calculation of AE'^^^ is needed, the computational cost is relatively modest, 
largely the inversion of an operator. 

One obtains a very natural interpretation of the second-order correction ( ^2[ ) by using this closure equation. 
Consider the generalized-Thomas- Fermi problem, Eq. ( |3^ ) , and suppose the external potential is slightly mod- 
ified by the quantity 5Vcxt{'<^)- One thus obtains a new solution of the GTF equation jt-qxf ~ '^gtf + Suq^f 
which would verify 



["GTf] + Kff ["-gtf] + (^Kxt = ^J■'GTF ■ (40) 



(57tf 
dn 

Subtracting Eq. (^) as before yields 

/ rfr^^(r,r')5nGTF(r') + J dr^(r, r') <SnGTF(r') + ^Kxt = A^ . (41) 

If we now choose the variation of the potential to be 

<5Fc.t(r) = Jdr' {S^£tot/Sn^)[nGTF]{r,r') h°'''{r') , (42) 

JriGTF + n"^"^ satisfies the same Eq. ( |39| ) as Sn. This means that, at this level of approximation, Sn is the sum 
of n°^'' and the displacement of charges Stiqtf screening n°^'^ in the GTF approximation. Indeed, the definition 
of a screened interaction V^c implies 

dr'^(r, r') Sn{r') = / dr'Vsc{r, r') n°'^{r') , (43) 



and therefore that the second-order correction, Eq. (32), can be written 



A£;(2) = i y drdr' n°"'=(r) V;c(r, r') n°""(r') . (44) 

Thus the second-order correction is simply the energy of interaction between the additional charge oscillations 
caused by the quantization, where the interaction is screened because, after all, the "other" electrons treated in 
GTF are around. Note that Vsc is the screened interaction within the finite sized system, not in the bulk, and 
so includes boundary effects |^ ; under certain conditions, the bulk screened potential may be used . More 
importantly, while the screened interaction here does include exchange-correlation at the GTF (macroscopic) 
level, the result ( ^ ) is a "direct-like" contribution while an "exchange-like" term is missed. This is related to the 
deficiencies of the LDA-like treatment of DFT here and presumably could be fixed through a local-spin-density 
functional approach. Such generalizations would be straight forward using the same arguments given here. 
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IV. CONTRIBUTION OF THE RESIDUAL INTERACTION TO PEAK SPACING DISTRIBUTIONS 



As an example of the utility of the Strutinsky method for adding quantization effects to a macroscopic result, 
we turn to considering the spacing between peaks in the conductance through a quantum dot in the Coulomb 
blockade regime. The contribution of residual interactions have been estimated for chaotic systems within a 
random matrix theory framework |28 4^-^. There it was found to be small, but not too far from the scale 



necessary to explain the failure of the constant interaction model. Here our ultimate aim is to evaluate the 
effect of the residual interaction in specific model systems which often are not in a regime where their quantum 
properties have fully converged to the purely statistical behavior found in random matrix theory. Systems tend 
not to be purely chaotic, and even when chaotic, still exhibit manifestations of short time dynamics in their 
eigenproperties. This can often lead to important deviations from statistical limiting behaviors. We therefore 
briefly sketch the relationship between the residual interaction and the Coulomb blockade peak spacings. 

The position of a conductance peak as a function of gate voltage is proportional to the change in the total 
energy of the system when an electron is added , 

^iN = E{N) - E{N - 1) , (45) 

and the conductance peak spacing is proportional to the discrete inverse compressibility 

Xn = /iw+i - ^J^N (46) 
= E{N + l) + E{N -1)-2E{N) . (47) 

For each of the ground state energies here we will insert the second-order Strutinsky approximation to the DFT 
energy. The first term, £'gtf7 is the ground state energy in the generalized-Thomas- Fermi approximation, and 
is essentially the charging energy of the dot. The first-order correction contains the single-particle quantization 
effects. In some sense these two terms together constitute the same level of approximation as the much used 
constant interaction model. In fact, more physics is included here since changes in the self-consistent confining 
potential |^ are explicitly contained in the Strutinsky approach , whereas due to the ad hoc nature of the 
constant interaction model, therein exists no information at all on the self-consistent potential. The second-order 
correction term, Ai?(^\ contains, then, the effects of the residual interaction. 



V. THE QUARTIC OSCILLATOR: A CASE STUDY 

Let us now illustrate the above approach with a particular example. For the sake of simplicity, we choose 
a one-dimensional model system consisting of N electrons in the confining potential Vext(2^) = a;^/2 with the 
interactions governed by the one-dimensional Poisson equation cP'Vi-at[n]{x) / dx^ = — 47re^n(a;). This is a simple 
limit of a three-dimensional problem: the system is assumed to be invariant in the tranvserse directions y and 
z so that the interactions are between planes of charge, but the medium is extremely inhomogeneous with the 
transverse mass taken to infinity so that only one-dimcnsional quantum mechanics is needed. Exchange and 
correlation effects are turned off; thus the interaction functional is 

/oo 
n{x')\x-x'\dx' . (48) 
-oo 

Note that use of the ID Poisson equation causes an interaction which grows with distance. Use of the subscript 
"int" in this section, rather than "coul" above, is meant to distinguish this case from the three-dimensional 
Coulomb interaction. We emphasize that our interest in this simple model system is only as an illustration for 
better understanding of the Strutinsky method. 

We vary the electron charge e to see how well the Strutinsky scheme works for different strengths of the 
interaction, e — 0.5, 1.0 and 1.5 in units where Ti — m — 1. The electron spin degeneracy is not considered here. 
First, we perform generalized-Thomas- Fermi and density- functional-theory calculations directly. Next, using the 
GTF results, we apply the Strutinsky techniques to find approximate DFT results. Finally, these approximate 
results are compared to the actual DFT values. Because of the neglect of exchange-correlation here, the GTF 
approach reduces to true Thomas-Fermi and the DFT approach is simply the coupled Schrodinger-Poisson 
equations. 
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A. Thomas-Fermi Numerical Calculations 



For one-dimensional systems 



and, thus, the kinetic energy term, Eq. (Q), can be written explicitly as 

2t2 ^oo 

TtfN = / n{xfdx. (50) 

The ground state density is obtained by solving the Thomas-Fermi equation [cf. Eq. (|^)] 



2m 



GTf(2;) + + Vint[nGTF]{x) = fJ-GTF (51) 



where we have used (STtp / Sn)[n\ = e(n) = {n'^h^ /2m)n{x)'^ . By differentiating twice and using the Poisson 
equation, one obtains the second order differential equation 

This can then be transformed into coupled first-order equations 

yi{x) = nl.rj,p{x) 
dyi{x) 



dx 

dy2{x) 2m 



y2{x) 

(W 7^^-6x2) (53) 



dx TT^h^ 

which can be conveniently solved. Because of the symmetry of the system, dn/dx = at the origin and one need 



only specify the density at the origin as an initial condition. One repeats solving Eq. (52) adjusting nQT^p(x = 0) 
on each iteration until the normalization condition N — J nQTp{x)dx is satisfied. Once the electron density 
noTvix) is found, the ground state energy is obtained from 

^^GTF = "^TF ["-GTf] + ^cxt ["-GTf] + ^int ["-GTf] (54) 

where 7tf is given in Eq ( ^2|) and 

f°° 1 

foxt["-GTF] = 2 / nGTFix)-x'^dx 
Jo ^ 

finthGTF] = 2 • - / nGTF(a;)Vlnt["-GTF](a;)rfa; ■ (55) 
^ Jo 

The electron densities ncTF for = 5, 10 and 20 with e = 1.0 are plotted in Fig. 1(a) and the effective 
potential, T4ff [tigtf] given by Eq. (||), in Fig. 1(b). All three cases show the same basic structure which can be 
simply understood as follows. Without the interaction (e = 0), the density would have one maximum at the 
origin since the external potential has a minimum at the center. Once the interaction is turned on, electrons 
repel each other and avoid the center, making two maxima in the density. Though not pictured, the larger the 
value of the electron charge e, the lower the central valley in the density, and the more the density maxima 
move away from the origin. As intuitively expected, the minimum points in the effective potential correspond 
to the maximum points of electron density, and increasing e increases rapidly the bimodal nature of the density. 
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B. Quantum Numerical Calculations 



The numerical calculation of the DFT energy requires self-consistently solving 

- 2^ + Vcs [n] {x)j (f)i {x) = ei4)i {x) 

N 



i(x)=^|0,(x)p (56) 



/oo 
n{x)\x-x \ dx 
-oo 

which are the coupled Schrodinger-Poisson equations. We start the self-consistent iterations with the Thomas- 
Fermi potential Vi3ff[nGTF]- At each iteration, we first diagonalize the Hamiltonian H = fy^ /2m + 14ff(r) 
expressed in the basis of Hq = fsP' /2+x^ /2. From the eigenvalues and eigenvectors of iJ [V^ff [n]] , we can construct 
the electron density and the corresponding effective potential. Self-consistency is evaluated by comparing the 
effective potentials V^'^ and before and after each iteration (or equivalently the densities n^'*^ and n'^®'^). 

Because of the well-known instability of the Poisson equation, one cannot simply use the output from one 
iteration, V^g''^, as the input to the next |46j. Instead, we feedback only part of the output 

^Old, Next Iteration ^ -^Old ^ ^^y^^.^ _ -^Old)^ Q < a < 1 (57) 

where a is initially set as 0.5. If the self-consistency is not improved, a is reduced by half and the iteration redone 
so that improvement is guaranteed for every iteration. We repeat this until the density reaches self-consistency, 

max |n'^™(x) - n'^^^{x)\ < 10"^ (58) 

We require self-consistency in the density rather than the potential because the overall magnitude of the density 
does not change significantly as N increases. 

Once the self-consistent density and effective potential are obtained, the Weyl part of the density, n^rp, as 
well as the chemical potential /xJ^ft can be calculated from Eq. (|l5|) ; the oscillating part of the density follows 
from n^prp ~ TiDFT — '^DFT- Finally, the self-consistent ground state energy is obtained using 



'?DFt["DFt] = flp[Kff[".DFT]] - 2 / nBFTix)Vcs[nBFT]{x)dx 

Jo 



(59) 



and the same expressions for fext and fint as in the Thomas- Fermi calculation, Eq. (pq). We have used the 
above relation for the kinetic energy instead of the definition since the eigenvalues are more stable than the 
eigenvectors in the numerical calculations. 

The quantum electron densities for N — 5, 10, and 20 are superposed in Fig. 1(a) for electron charge e — 1.0. 
One can see the quantum mechanical oscillations whose number of peaks corresponds to the electron number 
N. Note the decreasing oscillation amplitudes with increasing particle number, as well as the tunneling outside 
of the potential wall at the classical turning points. The effective potentials, superposed in Fig. 1(b), are 
indistinguishable from the corresponding Thomas- Fermi potentials. 



C. Strutinsky Energy Corrections 

In order to find the approximate ground state energy using the Strutinsky method, we start with the Thomas- 
Fermi density and potential, calculated above, and quantize in this potential by solving the Schrodinger equation 

(~imS^ + Kff["'GTF]^ (t>i — ^i4>i for the eigenvalues and eigenvectors. 
The first-order energy correction is given by Eq. (|3^); in our example, the expression for the Weyl part reduces 

to 

fi^[Kff[nGTF]] = 2^ j {fi^ + 2Fefr[nGTF](a:)) ^/i^ - 14ff[nGTF](a;) dx . (60) 
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The second-order energy correction, from Eq. ([3^), is 

/ h°'"'{x)\x-x'\Sn{x')dxdx' 

oo ^ — oo 

= ATie^ I n°''%x)\ I {x - x')5n{x')dx'\ dx . (61) 



Jo VJx ) 

The required input follows from Eq. (p9|). This equation can be simplified by noting, first, for the kinetic 
energy term 

[ncTF] (a;, x') - ^ [u^tf] Six - x') . (62) 
Second, for the term depending on Vcs note that 

-^[nG'TF\{x,x')5n{x')dx' = -2TTe^ I 5n{x')\x - x'\dx' = Vn,t[5n]{x) (63) 



implies 



d^Vi,,t[5n]{x) 



dx"^ 

Thus, by taking the second derivative with respect to x of Eq. (^), one obtains 

TT^h^ d^ 



47re^(5n(x) . (64) 



TO dx^ 



{noTFix) ■ {Sn{x) - n°"'=(x)) } - 4TTe^5n{x) = . (65) 



This equation can be converted into the coupled first-order equations 

yi{x) = nGTF(a;) • {5n{x) - rT^'^x)) 
dyi{x) 



dx 



V2{x) 

"'^'""^ '°''=(a;)) (66) 



dy2 jx) _ 4TOe^ / yi{x) ^ 
dx nh'^ \nGTF{x) 

which can be conveniently solved. 

With the energy correction terms calculated, the Strutinsky scheme allows us to approximate the quantum 
ground state energy using essentially classical Thomas- Fermi quantities. We plot AE = E'dft — ^'Gtf, AE — 
AE'-^\ and AE - AE'-^^ ~ AE^'^^ as functions of N to see the series convergence of the Strutinsky scheme 
in Figs. 2, 3 and 4 for e = 0.5,1.0 and 1.5 respectively. In the first two cases, the convergence seems to be 
working well: for e = 0.5, without correction terms the error is of order 0.1 and smooth, while the first-order 
correction term improves the accuracy to 0.001, and the second-order term to roughly 0.0001. For e = 1.0, 
without correction terms the error is about 0.1, with the first-order correction the error reduces to 0.01, and 
with the second-order to 0.0004. 

For e — 1.5, before comparing the order of magnitude of the different terms, it is useful to say a few words 
about the odd-even structure clearly visible in Fig. 4, and actually already apparent for e = 1.0 at low N 
in Fig. 3. The origin of this behavior is not related to spin (which has not been included) but can be readily 
understood by looking at the lower panel of Fig. 5, which shows the effective Thomas- Fermi potential K.ff [uctf] 
for N = 20, e — 1.5. Here, one sees that the latter has developed a barrier at the center of the well, the top of 
which lies very close to the chemical potential. For the quantum case, this means that the one particle levels 
below the Fermi energy are organized as quasi-doublets. This naturally leads to an odd-even effect since for 
N even (odd), the last occupied orbital has an energy consistently below (above) that suggested by the Weyl 
approximation. 

Moreover, because the last occupied orbitals are close to the chemical potential and so near the top of the 
barrier, it is clear that semiclassical approximations will be "at risk" here. This is clearly seen for instance in 
Figs. 6 and 7, which, for coupling e — 1.5 and = 19 and 20 particles, shows a comparison between the exact 
(5n(r) = nDFT(r) — ''^GTF(r) and its approximation obtained using Eq. (p9). The two curves are almost on top 
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of each other everywhere, except in the middle of the well — that is, near the maxima of the barrier. In addition, 
one observes that in that place, the approximation is worse for an odd than for an even number of particles. 
This can be explained by the fact that in the former case the last occupied orbital is symmetric and thus has 
a probability maximum at the origin, while in the latter case the last occupied orbital is antisymmetric and 
so goes to zero. As a consequence the errors in the Strutinsky approximation scheme also display an odd-even 
structure clearly seen in Fig. 5. If, however, one concentrates on the even case, where the effect of the central 
barrier is less important, without correction terms the error is about 0.2, with the first-order correction the 
error is reduced to 0.02, and further reduced to about 0.003 if one includes the second-order corrections. 



D. The Peak Spacings 



Turning now to the inverse compressibility xn introduced in Sect. |rV|, we observe the same trend as for 
the groundstate energy. The approximation x^^^"^ , calculated strictly by the Thomas- Fermi approximation, 
already gives an excellent relative precision. Moreover, for e = 0.5 and e — 1.0, each term in the Strutinsky 
development significantly reduces the error (i.e. an order of magnitude or more for N = 18). For these cases 
however the approximation is already much better than the mean spacing A if the first correction is included. 
We therefore show the data only for e = 1.5, for which the corrections are enhanced by the proximity of the 
chemical potential to the top of the inner barrier. In Fig. 8(a), the xn are shown. The discrete points represent 
the full quantum calculations and are taken as the reference points. It is seen that x^'^''"'' does not capture 
the odd-even double-well effect, but otherwise captures the essential peak spacing behavior. In Fig. 8(b), the 
relative errors are shown as a function of N . More specifically, the difference between the quantum xn and 
one, two, or three terms of the Strutinsky series — X^^''''^ X^^''^'' and XAr'"^^^\ respectively — is divided by A, 
the mean level separation. It is seen that the majority of the odd-even, double-well effect is captured by the 
first correction term. Reassuringly, even in this case for which the inner barrier degrades the quality of the 
semiclassical development, each inclusion of an extra term in the series reduces somewhat the relative errors of 
the Strutinsky method. Moreover the improvement due to the addition of the second order correction becomes 
more significant with increasing N . Thus we see that even in this more difficult case, the Strutinsky method 
gives a excellent scheme of approximation in the semiclassical limit. 



VI. CONCLUSION 



We have developed an approximate series expansion for studying the ground state of interacting systems 
using the idea of the Strutinsky shell correction method. We tested the validity of the Strutinsky scheme by 
numerical calculations of interacting electron systems in a one-dimensional, externally applied quartic potential. 
By varying the electron charge strength, we were able to confirm the stability of the method. It approximated 
extremely well the quantum mechanical DFT quantities using [semi] classical Thomas-Fermi data for three 
different charge strengths. One exceptional circumstance giving less reliable results was noted with respect 
to a barrier in VcR approaching the chemical potential. The calculations show a tendency for the series to 
converge as electron number increases. This is consistent with expectations for systems with large numbers of 
particles because of the increasing reliability of the Thomas-Fermi calculations as the system goes deeper into 
the semiclassical regime. 

The method discussed in this paper could serve two purposes. On one hand, it gave us an efficient way 
to proceed with numerical calculations, and it is conceivable that this approach could be of some help for 
larger scale, realistic DFT calculations [|4 35 3^. On the other hand, it also provides some physical insight by 
decomposing the total DFT energy into various contributions, each of them receiving an intuitive interpretation. 

In the context of quantum dots, for instance, one of its simple but presumably useful applications could be to 
justify, and make more precise, the constant interaction model. Indeed, this latter is usually presented as an ad- 
hoc model motivated essentially by its simplicity. Here, up to some reinterpretation of what is the capacitance 
of the dot, we see that the constant interaction model can be understood as the first-order approximation in 
a Strutinsky development of a DFT calculation. One obtains, in addition, that the various parameters of the 
model (classical energy and potential governing the motion of the independent particles) are specified, and, in 
principle, can be computed explicitly. This makes it possible, for instance, to study the sensitivity of the dot's 
energy to the variation of an external parameter |^5|] . 
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In the same way, the second-order correction term gives insight into the "residual" interactions between 
electrons. In the context of DFT it gives some basis to the fact that electrons in quantum dots behave as Landau 
quasi-particles interacting weakly through a screened Coulomb interaction. It moreover explicitly specifies how 
the screening process is affected by the confinement of the electrons, which might be relevant in the limit where 
the screening length is not much smaller than the size of the dots. 

With the failure of the theoretical predictions of the single-particle constant-interaction model to explain the 
experimentally observed conductance peak spacing distributions |17-2^], our original interest was to investigate 



what physical factors are missing and to understand better the statistical behaviors of quantum dots. We expect 
that the approach of defining the many-body system, explicitly including electron-electron interactions within 
density functional theory, to bring us closer to a resolution of the failure of the constant interaction model in 
this context. We leave application to more realistic systems for future work. 
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APPENDIX A: 

Consider a Hamiltonian H — p^/2to + V{y) and its perturbed Hamiltonian H' = H + 6V{r). We denote the 
eigenvalues and eigenvectors as and 0i(r) [e^ and (/''(r)] for H [H'). 
To second order in 5V , the perturbed eigenvalues are 



^, = e., + e^+^ (Al) 



where 



^(2) ^ mm'i^3)\ 

^ ^-^ f — f 

... tj t 



Y^ mwm'' ^ (A2) 



assuming non-degenerate eigenstates. 

Similarly, taking H' as the original Hamiltonian and H = H' — 5V as the perturted Hamiltonian, one can 
write 

= (A3) 

where 



e. 



P=6f'+ 0(51/3). (A4) 



Subtracting Eq. ( [A3| ) from Eq. ( Al ) , one obtains to second order in SV 



= \mm'l^^) + mm<p'^))■ (A5) 



Summing over i, one obtains 

1 ^ 



1=1 



= \ I SV{r){Ti{r) + n'{r))dr . (A6) 
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FIG. 1. The electron densities (upper panel) and effective potentials (lower panel) for interacting particles in a quartic 
potential. The results for both the Thomas-Fermi (dashed) and coupled Schrodinger-Poisson (solid) approximations are 
given. The electron charge is e = 1.0, and the electron number is iV = 5, 10, and 20 from bottom to top in upper panel, 
top to bottom in lower panel. V^ff [noTp] and V"efi[nDFT] coincide so well that one cannot distinguish the differences here. 
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FIG. 2. Convergence of approximations to the quantum ground state energy for electron charge e — 0.5. The curves 
are, from top to bottom, the error in the Strutinsky energy correction approach taken at zeroth, first, and second order: 
specifically AE = Eoft - Sqtf (dotted), AE - AE'-^^ (dashed), and AE - AB^^^ - AE^^^ (solid). The convergence in 
this case is excellent. 
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FIG. 3. Convergence of approximations to the quantum ground state energy for electron charge e = 1.0. The curves 
are, from top to bottom, AE = Euft - Egtf (dotted), AE - AE^^^ (dashed), and AE - AE^^^ - AE^^'^ (solid). 
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FIG. 4. Convergence of approximations to the quantum ground state energy for electron charge e = 1.5. The curves 
are, from top to bottom, AE = Eoft - Eqtf (dotted), AE - AE'^^'> (dashed), and AE - AE'^^'> - AE^'^'^ (solid). 
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FIG. 5. The electron density (upper panel) and effective potential (lower panel) for N = 20 particles and e = 1.5. In 
the upper panel, results for both the Thomas-Fermi (dashed) and quantum (solid) cases are given. In the lower panel, 
the dashed horizontal line is the position of the chemical potential /^tf . 
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FIG. 8. The discrete inverse compressibility as a function of electron number A'^. The open circles are xn using the 
quantum ground state energy (exact). Results using three approximate ground state energies are shown: (1) dotted: 
Thomas-Fermi, x^^*''' using -Egtp, (2) dashed: first-order Strutinsky, x^^*^' using -Egtf + AE'-^\ and (3) solid: 
second-order Strutinsky, x^'^'^' using _Egtf + AiS*-^^ -I- Ai?''^^ The lower panel shows the relative errors, (xiv — x^^)/A, 
of the same three approximations. 
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